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Abstract: We apply the positive-weight Monte Carlo method of Nason for simulating 
QCD processes accurate to Next-To-Leading Order to the case of e + e~ annihilation to 
hadrons. The method entails the generation of the hardest gluon emission first and then 
subsequently adding a 'truncated' shower before the emission. We have interfaced our 
result to the Herwig++ shower Monte Carlo program and obtained better results than 
those obtained with Herwig++ at leading order with a matrix element correction. 
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1. Introduction 

Matching next-to-leading order (NLO) calculations to shower Monte Carlo (SMC) models 
is highly non-trivial. One method (MC@NLO) has been proposed in [1] and has been 
successfully applied to several processes [2, 3] in connection with the FORTRAN HERWIG 
parton shower. One drawback of the method is the generation of negative weighted events 
which are unphysical. 

An alternative method proposed in [4] overcomes the negative weight problem. It has 
successfully been applied to Z pair hadroproduction [5] and involves the generation of the 
hardest radiation independently of the SMC model used to generate the parton shower. To 
preserve the soft radiation distribution, the addition of a 'truncated shower' of soft coherent 
radiation before the hardest emission is necessary. In this report, we will be implementing 
this method for electron-positron annihilation to hadrons at a centre-of-mass energy of 91.2 
GeV. The SMC we will be employing is Herwig++ [6] which will perform the rest of the 
showering, hadronization and decays. 

In Section 2, we discuss the generation of the hardest emission according to the ma- 
trix element via a modified Sudakov form factor. In Section 3, we go on to describe the 
generation of a simplified form of the 'truncated' shower. We determine the probability for 
the emission of at most one gluon before the hardest emission and re-shuffle the momenta 
of the outgoing particles accordingly. This is a reasonable first approximation since the 
probability of emitting an extra gluon is small. 

In Section 4, we present our results compared with the measured data and those 
obtained using Herwig++ with a matrix element correction. Finally in Section 5, we sum- 
marize our conclusions. 
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2. Hardest emission generation 



2.1 Hardest emission cross section 



The order-a s differential cross section for the process e + e 
masses, can be written as 



qqg, neglecting the quark 



R{x,y) = a W(x,y) = a 



2a, 



x 2 + y 2 



3vr (l-x){l-y) 



(2.1) 



where do is the Born cross section, W = R/<Jq, x and y are the energy fractions of the 
quark and antiquark and 2 — x — y is the energy fraction of the gluon. R has collinear 
singularities when the quark or antiquark is aligned with the gluon. When combined with 
the virtual corrections, these singularities can be integrated to give the well known total 
cross section to order-a s , 



a, 
1 + — 

7T 



pnlo = &o 

From [4], we can write the cross section for the hardest gluon emission event as 



da = ^2,B(v)d$ v 



(2.2) 



(2.3) 



where B(v) is the Born cross section and v is the Born variable, which in this case is the 
angle 6 between the e + e~ beam axis and the qq axis, r represents the radiation variables 
(x and y for our specific case), d& v and d& r are the Born and real emission phase spaces 
respectively. 

A^ lo (pt) is the modified Sudakov form factor for the hardest emission with transverse 
momentum px, as indicated by the Heaviside function in the exponent of (2.4), 



\NLO 
^R 



(Pi 



cxp 



d*r^^-@(k T (v,r)- PT ) 
B(v) 



Furthermore, 



B(v) = B{v) + V(v) + J (R(v, r) - C(v, r))d<S> r . 



(2.4) 



(2.5) 



B(v) is the sum of the Born, B(v), virtual, V(v) and real, R(v,r) terms, (with some 
counter-terms, C(v,r)). It overcomes the problem of negative weights since in the region 
where B(v) is negative, the NLO negative terms must have overcome the Born term and 
hence perturbation theory must have failed. Now explicitly for e + e~ annihilation, 

,2 



a NLO 



(pt) 



exp 



dx dy- 



2a, 



x 2 +y' 



-Q(k T (x,y)) - pt) 



3vr (l-x)(l-y) 

where x and y are the energy fractions of the quark and antiquark and we define 



kr(x,y) = 



(1 - x)(l - y)(x + y - 1) 



max(x, y) 2 



(2.6) 



(2.7) 



and s as the square of the center-of-mass energy, kx is the transverse momentum of the 
hardest emitted gluon relative to the splitting axis, as illustrated in Figure 1 below. 
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Figure 1: Transverse momentum, kr- 
2.2 Generation of radiation variables, x and y 

From (2.1) and (2.3), it can be deduced that the radiation variables are to be generated 
according to the probability distribution 

A w {k T )W{x,y)dx dy 



(2.8) 



where in the particular case of e + e annihilation, A w (kx) and W(x,y) are given in (2.6) 
and (2.1) respectively. However for ease of integration, we will use a function 



U(x,y) = 



2a, 



3vr (l- x )(l-y) 



> W(x,y) 



(2.9) 



in place of W(x,y). As we shall see later, the true distribution is recovered using the veto 
technique. The variables x and y are then generated according to 



A u (k T )U(x,y) = U(x,y)exp 
This is outlined in the following steps: 

i) Set p max = k T max- 



J U(x,y)G(k T (x,y) - p T )dxdy 



(2.10) 



ii) For a random number, n between and 1, solve the equation below for px 

A U ( PT ) 



n = 



A u (p max ) ' 

iii) Generate the variables x and y according to the distribution 

U(x,y)5(k T (x,y) - p T ) ■ 



(2.11) 



(2.12) 



iv) Accept the generated value of pt with probability W/U. If the event is rejected set 
Pmax = Pt and go to step ii). 
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Now for e + e annihilation, 

•1 /•! 



U(x,y)Q(k T (x,y) - p T )dxdy = dx dyU(x,y)Q(k T (x,y) - p T ) ■ (2.13) 

J JO Jl-x 

In the region where x > y, let's define the dimensionless variable, k as 

K = 4 = ^-x){l-y){x + y-l) 
s x 2 

There are 2 solutions for y for each value of x and k, i.e y\ = y and yi = 2 — x — y. 



-2 + 3x - x 2 =F Vx 2 - 2x 3 + x 4 - 4kx 2 + 4kx 3 ,„ 
yi ' 2 = 2^1) • ^ 

Exchanging the y variable for n in (2.13) we find 

U(x,y)Q(hr(x,y) - p T )dxdy = J dx J dn — s — — —U(x,k) 



dx / 

- min K 



dn 



3ir (1 — x + p(n, x))p(n, x) 

(2.16) 



where 

p(n, x) = v / (l-2;)(l-4«;-x) . (2.17) 

Note that the argument of a s is the scale fop 2 = ks with Aqcd = 200 MeV. Equation 
(2.16) applies to the region of phase space where x > y. For the region x < y, x and y are 
exchanged in the equation. 
The x integration can be performed to yield 



U(x,y)@(k T {x,y) - p T )dxdy = J dn 



2a s (ns) 
3vr 



ln(l — x) 



21n(VT^x + VI - 4k - x) 



(2.18) 



For k < 0.08333, the region of phase space where x > y and the two solutions for y for 
given values of x and k are illustrated in Figure 2 below. The two solutions lie on either 
side of the dashed line y = 1 — ^x and are equal when x = x max = 1 — 4k which lies on 
the dotted line. At k = 0.08333, the branches meet along the line y = x and there is only 
one solution for y in the region (the lower branch). So for k > 0.08333, only one y solution 
exists. In addition there are no y solutions for k > 0.09. This is illustrated in Figure 3. 
Also note that for x < x u , there is only one solution for y. 

In the region where there are two solutions, the integral in (2.18) is performed along 
both branches independently and summed. For the upper branch, x runs from x u to 
^max = 1 — 4k while for the lower branch, x runs from x l to x max where if we define 
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x a = 39k - 1 + k 3 + 15k 2 , 
Xb = 33k 2 + 3k — 3k 3 



x c — \l x a ^ -\- Xb^ 

, -1 f x b 
Xd = tan — 

V X n 



X& — 



x f 



-i2 Xc COS (y) ' 

(-1-k 2 -10k) cos(^) 



12x1 



K + 5 




X h = :-ui ("■-'■ ) j r;< H ■ ^— ) . (2..19) 

we can write x u and x l as; 

X — Xp ~\~ Xy ~h X^f ~h X/^, , 

x' = x e + Xf + Xg - Xh] (2.20) 

In the region where there is only one solution for y, x runs from x l to x u . 
The k integration can then be performed numerically. Having performed the integra- 
tion, values for k and hence Ut are then generated according to steps 1 and 2 in Section 
2.2. In step 3 of Section 2.2, the variables x and y are to be distributed according to 
U(x,y)5(kT{x,y) —pr)- This is the subject of the next section. 

2.3 Distributing x and y according to U(x, y) 

To generate x and y values with a distribution proportional to U(x,y)5(kT(x,y) — pr), we 
can use the 5-function to eliminate the y variable by computing 

U(x,y) 



D(x) = / 5(k T -p T )U(x,y) 



dy 



(2.21) 



y=y 



where y is such that kT{x,y) = pt- Note that is the same for both y solutions. We 
then generate x values with a probability distribution proportional to D with hit-and-miss 
techniques as described below. All events generated have uniform weights. 

i) Randomly sample x, N x times (we used N x = 10 5 ) in the range [a; m i n : x max ] for the 
selected value of k. 



ii) For each value of x, evaluate D = D(x, y\) + D(x, y-z) if there are two solutions for the 
selected k and D = D(x,y2) if there is only one solution. Also, if k < 0.08333 and 
x < x u (see Figure 2), there is only one y solution so evaluate D = D(x,y2)- 
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iii) Find the maximum value D mSuX of D for the selected value of k from the set of N x 
points that have been sampled. 

iv) Next, select a value for x in the allowed range and evaluate D. 

v) If D > rl) max (where r is a random number between and 1), accept the event, oth- 
erwise go to iv) and generate a new value for x. 

vi) If for the chosen value of x, there are two solutions for y, select a value for y in the 
ratio D(x,yi) : D(x,y 2 ). 

vii) Compare U(x,y) with the true matrix element, W(x,y). If the event fails this veto, 
set K max = k and regenerate a new k value as discussed in Section 2.2. 

NB: For the region y > x, exchange x and y in the above discussion. In this way, the 
smooth phase space distribution in Figure 4 below was obtained for the hardest emission 
events. The plot show 2,500 of these events. 
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Figure 4: Phase space and distribution of hardest emissions. 
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3. Adding the truncated shower 



As mentioned in Section 1, a 'truncated shower' would need to be added before the hardest 
emission to simulate the soft radiation distribution. Due to angular ordering, the 'trun- 
cated' radiation is emitted at a wider angle than the angle of the hardest emission as well as 
a lower pr- This means the 'truncated' radiation does not appreciably degrade the energy 
entering the hardest emission and justifies our decision to generate the hardest emission 
first. 

Below is an outline of how the 'truncated shower' was generated. We will consider 
the case in which at most one extra gluon is emitted by the quark or antiquark before 
the hardest emission. The outline closely follows the Herwig++ parton shower evolution 
method described in [6, 7] where the evolution variables z, the momentum fractions, and 
q, the evolution scale, determine the kinematics of the shower. 

i) Having generated the pr of the hardest emission and the energy fractions x and y, 
calculate the momentum fraction z and 1 — z of the partons involved in the hardest 
emission. We will assume henceforth that x > y and that y is the energy fraction of 
the quark, i.e. the quark is involved in the hardest emission. Then 



where as defined above z and 1 — z are the momentum fractions of the quark and gluon 
respectively. 

ii) Next generate the momentum fraction z% of the 'truncated' radiation according to the 
splitting function P qq = within the range 



where q~i is the initial evolution scale, i.e. yfs = 91.2 GeV, /i=max(m a , Q g ), m a is the 
mass of the quark and Q g is a cutoff introduced to regularize soft gluon singularities 
in the splitting functions. In this report, a Q g value of 0.75 GeV (and hence [i = 0.75 
GeV) was used. z% is the momentum fraction of the quark after emitting the 'trun- 
cated' gluon with momentum fraction 1 — Zf. 

iii) Having generated a value for z t , determine the scale q\ of the hardest emission from 




(3.2) 




(3.3) 



where z = Zt- 
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iv) Starting from an initial scale q~i, the probability of there being an emission next at the 
scale q is given by 



S(q~i,q) 



A(? c ,g) 



where 



j dz^P qq 9(0<p T <p T ) 



(3.4) 



(3.5) 



g c is the lower cutoff of the parton shower which was set to 0.4 GeV in this report, a s 
is the running coupling constant evaluated at z(l — z)q, P qq is the q — ► qg splitting 
function and kx is the transverse momentum of the hardest emission. The Heaviside 
function ensures that the transverse momentum, p T of the truncated emission is real 
and is less than pr- To evaluate the integral in (3.5), we overestimate the integrands 
and apply vetoes with weights as described in [6]. With r a random number between 
and 1, we then solve the equation 



S(qi,q) 



(3.6) 



for q. If q > q^, the event has a 'truncated' emission. If q < q\ , there is no 'truncated' 
emission and the event is showered from the scale of the hardest emission. 



v) If there is a 'truncated' emission, the next step is to determine the transverse momen- 
tum p T of the emission. This is given by [6] 

A = yJ(l-z t y^-^)-z t Q g 2 . (3.7) 
If p T 2 < or p T > pt go to ii) . 



vi) We now have values for zt, the momentum fraction of the quark after the first emission, 
p T , the transverse momentum of the first emission, z, the momentum fraction of the 
hardest emission and pr, the transverse momentum of the hardest emission. This 
is illustrated in Figure 5. We can then reconstruct the momenta of the partons as 
described in [6]. The orientation of the quark, antiquark and hardest emission with 
respect to the beam axis is determined as explained there for the hard matrix element 
correction. 



In principle this procedure could be iterated to generate a multi-gluon truncated 
shower. However, for the present, we consider only the effect of at most one extra gluon 
emission. As discussed in [4], the initial showering scale of the hardest gluon is set equal 
to the showering scale of the quark or antiquark closest in angle to it. This gives the right 
amount of soft radiation colour connected to the gluon line. 
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Figure 5: Adding the 'truncated' emission. (E = 45.6 GeV) 

4. Results and data comparisons 

One million events were generated as described above and then interfaced with the SMC, 
Herwig++. 13% of the events acquired an extra 'truncated' gluon. A px veto was imposed 
on the subsequent shower starting from the hardest emission to the hadronization scale 
which was tuned to 0.4 GeV. Table 5 and Figure 6 show respectively the average multiplic- 
ities of a wide range of hadron species and the charged particle multiplicity distribution. 
The subsequent figures are plots of comparisons with event shape distributions from the 
DELPHI experiment at LEP [8]. 

The upper panel below the main histograms shows the ratio Mi D D% (where Mj and 
Di stand for Monte Carlo result and data value respectively) compared with the relative 
experimental error (green). The lower panel shows the relative contribution to the x 2 °f 
each observable. As in [6], the x 2 contributions allow for a 5% uncertainty in the predictions 
if the data are more accurate than this. Finally, in Table 2 we show a list of x 2 values 
for all observables that were studied during the analysis. The results were generated using 
Herwig++ 2.0 [9] which includes some improvements in the simulation of the shower and 
hadronization as described in [9] leading to some changes in the x 2 f° r specific observables 
relative to [6], although in most cases the changes are small. 

5. Conclusions 

We have successfully implemented the Nason method of generating the hardest emission 
first and subsequently adding a truncated shower for e + e~ annihilation into hadrons. 

We have tested the method against data from e + e~ colliders and for almost all ex- 
amined observables, the simulation of the data is improved with respect to Herwig++. In 
particular the Nason method seems to fit the data better in the soft regions of phase space. 
The poorer fits obtained for variables such as the thrust minor which vanish in the three-jet 
limit (and in general for planar events) may be attributable to the lack of multiple emission 
in the truncated shower. 
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The fits are summarized in the \ 2 values shown in Table 2. In Table 2 we also present 
the x 2 values for the observables obtained without the implementation of the 'truncated' 
shower. We see there is a general improvement in the fits associated with the addition of 
the 'truncated' emissions. 

Future work in this area will extend this method to processes with initial state radiation 
in hadron-hadron collisions aiming at simulating Tevatron and LHC events. 
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0.097 ± 0.007 
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P 
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0.0471 ± 0.0046 
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0.0262 ± 0.001 


0.032* 


0.034* 
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0.0058 ± 0.001 
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01 8zL 


D± 


A,D,0 


0.184 ± 0.018 


0.252* 


0.282* 


L>*(2010) ± 


A,D,0 


0.182 ± 0.009 


0.173 


0.185 


#° 


A,D,0 


0.473 ± 0.026 


0.497 


0.497 


Df 


A,0 


0.129 ± 0.013 


0.130 


0.150 


D? 





0.096 ± 0.046 


0.06 


0.044 


J/* 


A,D,L,0 


0.00544 ± 0.00029 


0.00349* 


0.0033* 


A+ 


D,0 


0.077 ± 0.016 


0.0215* 


0.0303* 


*'(3685) 


D,L,0 


0.00229 ± 0.00041 


0.00165 


0.00154 



Table 1: Multiplicities per event at 91.2 GeV. We show results from Herwig++ with a matrix 
element correction (Herwig++ ME) and the new method, Nason@NLO . Experiments are Aleph(A), 
Delphi(D), L3(L), Opal(O), Mk2(M) and SLD(S). The * indicates a prediction that differs from 
the measured value by more than three standard deviations. Data are combined and updated from 
a variety of sources, see ref [10]. 
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Figure 6: The distribution of the charged particle multiplicity. Data from the OPAL experiment 
at LEP [11]. 
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Figure 8: Oblateness, Sphericity, Aplanarity and Planarity distributions. Data from the DELPHI 
experiment at LEP [8]. 
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Figure 9: C Parameter and D Parameter distributions and the high, Mhi g h and low, M\ ow hemi- 
sphere masses. Data from the DELPHI experiment at LEP [8]. 
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Figure 10: The wide and narrow jet broadening measures -B max and B lri \ n and the difference in 
hemisphere masses, M^ie- Data from the DELPHI experiment at LEP [8]. 
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Table 2: x 2 /bin for all observables we studied. 
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